The governing equations used in the present LES are three-dimensional (3D) filtered compressible Navier–Stokes equations, including the following conservation equations for mass, momentum, energy, and species:

The density-based compressible solver astroFoam,^[31,35] which is developed based on the basic rhoCentralFoam solver^[38] distributed with OpenFOAM v2.3.0, is used to solve the above equations. Comparing with the original rhoCentralFoam solver, two aspects are modified to create the astroFoam solver. First, multiple species transport and multi-component diffusion (Eq. (4)) are added into the astroFoam solver to model more realistic mixing and reacting flows. Second, the energy equation (Eq. (3)) is solved for senible enthalpy in astroFoam instead of total energy in rhoCentralFoam to include the chemical reaction and species transport terms more easily. In addition, the convection–diffusion equation is solved by semi-discrete KT scheme^[39] for capturing shock and solving turbulence in astroFoam. A normalized variable diagram (NVD) scheme^[40,41] is used to reconstruct the primitive values at faces to obtain the second order accuracy. Time integration is carried out by the Crank–Nicholson scheme,^[41] which is second order accuracy in time.

Figure 1(a) shows the computational domain used in the current LES, which mainly consists of a box of size 50 mm × 100 mm × 50 mm in x, y, and z Cartesian coordinate directions. Previous studies^[6,7,11,16] indicated that the sound waves originating from the down stream will propagate upstream to change the initial shear layer structures at the nozzle exit, which will influence the development of the jet shear layer in the down stream further. However, the priori knowledge of nozzle exit conditions is usually difficult to obtain in many practical applications. In the present study, a convergent nozzle geometry is included in the upstream part of the computational domain. The nitrogen jet (with total pressure P₀ and total temperature T₀) is injected into the quiescent air (with static pressure P_∞ and static temperature T_∞) from the convergent nozzle of 20.0 mm in height. The entrance and exit diameters of the nozzle are d = 8.0 mm and D = 2.0 mm, respectively.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Simulation setup showing (a) computational model and (b)–(d) computational grid: (b) side view at z/D = 0, (c) top view at y/D = 0, (d) close-up view at nozzle exit.

As indicated by the previous studies,^{[27–31,42–46]} the spatial resolution in LES of supersonic jets needs to be rather high. Figures 1(b)–1(d) show the hexahedral, block-structured grid used in the present LES. The grid contains a refinement region of high resolution, covering the jet core and the jet shear layers. In the refinement region, two levels of grid resolutions, i.e., fine and coarse, are used to characterize the grid dependence of the solution. The fine and coarse meshes contain about 27.3 and 13.0 million cells, respectively, and their resolutions are similar to those used in previous LES of supersonic jets,^[27,43–46] which are summarized in Table 1. On the other hand, coarse cell sizes with a resolution of 1.0 mm in the far field and 0.5 mm at the outflow boundaries are used to avoid wave reflections by introducing additional dissipation.

Table 1.

Table 1.

Table 1. Grid resolution comparison in the near field of supersonic jets. . Grid D/mm Δrmin Δrmax Δymin Δymax Re Total/106 Fine 2.0 D/200 D/52 D/67 D/25 ∼105 27.3 Coarse 2.0 D/200 D/25 D/50 D/20 ∼105 13.0 Vuorinen et al.[27] 1.4 D/70 D/50 D/35 D/25 ∼105 12.0 Gorle et al.[43] 2.0 D/100 D/50 D/100 D/25 ∼105 17.4 Dauptain et al.[44] 25.4 D/35 D/30 D/35 D/30 ∼106 22.0 Rana et al.[46] 4.0 D/33 D/33 D/33 D/33 ∼104 9.2

Table 1. Grid resolution comparison in the near field of supersonic jets. .

The quiescent air is the mixture of 76.699% nitrogen and 23.301% oxygen in mass fraction. Initially the temperature, pressure, density, and velocity of the ambient air are all uniform, i.e., T_∞ = 300 K, P_∞ = 101325 Pa, ρ_∞ = 1.17 kg/m³, and U_∞ = 0. Four different simulations with NPRs of 5.60, 7.47, 9.34, and 11.21 are carried out to examine the shock behavior in underexpanded jets. The flow conditions at the nozzle exit are close to sonic conditions for different jets, and are summarized in Table 2 in detail. The stagnation condition for temperature and pressure is employed at the high-pressure nozzle inlet, while the velocity is treated with a zero-gradient condition. All walls are treated as no-slip, adiabatic condition. At the top of the computational domain and on the four free surfaces, an open boundary condition is used, i.e., all flow parameters are treated as the zero-gradient for outflow and fixed ambient values in the case of backflows.

The time step Δt in the LES of supersonic flows is generally small. Kawai and Lele^[45] used a time step satisfying Δt · U_∞/D = 9.6 × 10⁻⁴ in the LES of a sonic jet in supersonic cross flow. Liu et al.^[42] employed a maximum Courant–Friedrichs–Lewy (CFL) number of 0.4 while Dauptain et al.^[44] used 0.7 in their LES modeling of supersonic underexpanded jets. The time steps used by Hamzehloo and Aleiferis^[30] in the LES of methane underexpanded jets are both between 2 × 10⁻⁸ s and 5 × 10⁻⁸ s. The computational time step used in the present LES is similar to those used in previous LES studies,^{[30,42,44,45]} and takes a value approaching Δt = 1.37 × 10⁻⁸ s, which gives a Δt · a/D of 2.42 × 10⁻³ and is limited by the maximum CFL number of 0.6. The flow-through time (FTT) for jets washing out the computational domain in the streamwise direction is around 0.5 ms = 200t₀, with t₀ = 2.5 × 10⁻⁶ s = 2.5 μs. The total simulation duration is thus set to be 4FTT = 2.0 ms = 800t₀. The instantaneous results are saved every 2t₀, and turbulent statistics are collected for the last three flow-through times (3FTT, 200t₀ to 800t₀).

Table 2.

Table 2.

Table 2. Flow parameters and simulation conditions. . Property Symbol Case 1 Case 2 Case 3 Case 4 Units Mach number at nozzle exit M1 1.0 1.0 1.0 1.0 – Static pressure at nozzle exit P1 0.3 0.4 0.5 0.6 MPa Stagnation pressure P0 0.57 0.76 0.95 1.14 MPa Static temperature at nozzle exit T1 300.0 300.0 300.0 300.0 K Reynolds number at nozzle exit Re1 1.36 1.82 2.27 2.73 105 Static pressure ratio P1/P∞ 2.96 3.95 4.93 5.92 – Nozzle pressure ratio (NPR) P0/P∞ 5.60 7.47 9.34 11.21 –

Table 2. Flow parameters and simulation conditions. .

In order to establish the fidelity of the LES results, a grid convergence study has been conducted. The LES of underexpanded jet at NPR = 5.60 is carried out by using the two different grids, i.e., the coarse and fine meshes. It is important to capture the quasi-periodic shock structure as its interaction with the jet shear layer is responsible for the shock screech phenomenon. Figure 2(a) shows the mean pressure profiles along the jet centerline. As can be seen, the agreement over the first four shock cells is very good for the two meshes. The turbulence properties of the jets also need taking great concern. The radical profiles of the mean normal velocity, U_x, are compared for the two meshes in Fig. 2(b). The overall trends are similar for the two profiles, but the coarse mesh seems to, as expected, restrain the development of turbulence because of the larger grid dissipation, which is indicated by a relatively narrow jet boundary. Therefore, the fine mesh is mainly used for the results discussed subsequently, while the results obtained by the coarse mesh are only presented in Subsection 4.4.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Grid independence analysis for NPR = 5.60, showing (a) axial pressure profiles and (b) radical profiles of mean normal velocity at y/D = 10.

Figures 3(a) and 3(c) show the time-averaged density gradients for NPR = 5.60 and 9.34 obtained by LES, respectively. As can be seen, the typical wave structures in the near field of a highly underexpanded jet, including the Mach disk, intercepting shock, triple point, reflected shock, and slip lines that have been confirmed by previous experimental and numerical studies,^{[8,17,22–24,27–34]} are all well captured by the LES. The predicted shock structures are also in good agreement with those in the schlieren photographs shown in Figs. 3(b) and 3(d) under the same NPR. In particular, the locations of the Mach disk, the first and the second reflected shocks on the jet boundary (marked by lines) predicted by the LES are reasonably comparable to the measurements.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Comparisons of time-averaged near-field shock structure between LES and experimental data for NPR of (a), (b) 5.60 and (c), (d) 9.34.[35]

Figure 4 shows the iso-surfaces of density gradient for NPR = 5.60, which presents the 3D formation and development of the near-field shocks. The jet boundary, the intercepting shock, and the reflected shock are circular in cross-section. The slip line surface (or shear layer) initially expands the down stream of the Mach disk, then rapidly contracts as observed experimentally by Mitchell et al.^[23] The flow is supersonic (Ma ∼ 3.0) and possesses a low temperature (T ∼ 100 K) and high velocity (U_y ∼ 700 m/s≈ 2U_e, U_e is the nozzle exit velocity) in the vicinity of the Mach disk because of rapid expansion, and becomes subsonic (Ma ∼ 0.2) immediately after the Mach disk. Then the flow is accelerated to sonic conditions through the circular duct bounded by the shear layer.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Iso-surfaces of density gradient of 3D near-field shock structures for NPR = 5.60: (a) solid view, (b) perspective view.

Figure 5(a) shows quantitative comparisons of time-averaged density along the jet centerline between the LES results for NPR = 5.60 and the measured data by Panda and Seasholtz^[17] at a similar NPR of 5.74. Note that ρ_j is the fully expanded jet density. Good agreement with the experimental data is observed for the first two shocks. However, roughly from the third shock, the LES results start to deviate from the measurements. This deviation may be attributed to the difference between ρ_j/ρ_∞ values used in the LES and the experimental values. On the other hand, the radial profiles of density near the Mach disk obtained by the LES are compared with the measurements in Fig. 5(b). The LES results agree well with the measured data at the first two streamwise positions, i.e., y/D = 0.9 and 1.2. At y/D = 1.5, the predicted density peaks at around x/D = 0.2, which is smaller than the measured value of x/D = 0.25. This observation suggests that the present LES presents a narrower Mach disk than the experiment. Similarly, the LES studies performed by Vuorinen et al.^[27] and Hamzehloo and Aleiferis^[47] reproduced the width of the Mach disk smaller than the measurements. The contributing factors, as indicated by Franquet et al.^[8] and explained by Hamzehloo and Aleiferis,^[47] can be the differences in nozzle geometry design and simulation setup (e.g., the turbulence levels at the nozzle exit) between the LES modeling and the experiment.

In addition, several previous numerical studies^[48,49] have indicated the existence of a recirculation zone formed immediately downstream of the Mach disk in an underexpanded jet. However, the recent quantitative velocity measurements for a highly underexpanded jet at an NPR of 4.2 by Mitchell et al.^[23] showed that there is no experimental evidence for such a phenomenon. Figures 6(a) and 6(b) show the time-averaged profiles of pressure and streamwise velocities obtained by the present LES. As can be seen, the velocities for different jets indeed obtain the minimum values after the Mach disk, but there is no upstream velocity component (negative value) in the mean. This finding agrees reasonably with the prior experimental data,^[22,23] and also indicates that the present LES produces no recirculation region behind the Mach disk, which is in accordance with the recent LES study by Hamzehloo and Aleiferis.^[30]

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Time-averaged density comparisons between LES and the measurements by Panda and Seasholtz:[17] (a) axial profile, (b) radial profiles.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Time-averaged flow properties along jet centerline: (a) non-dimensional pressure, (b) streamwise velocity.

Figure 7 shows the time-averaged as well as instantaneous contours of density gradient on the midline plane at an NPR of 5.60 obtained by the LES. As also visualized through using schlieren photographs by Panda,^[15] the instantaneous flow fields are significantly different from the time-averaged ones, and undergo considerable distortion and show some oscillations. The distortion and oscillation progressively increase for the shocks formed further in the down stream from the nozzle exit. Especially at 5.0 < y/D < 7.5, the intense oscillation occurs, where the jet boundary presents a left and right wavy motion associated with the deformation of the rhombus shape of the third shock cell. These observations are two-dimensional (2D) impressions of a 3D helical motion of the jet where the jet core precesses about the jet centerline. After around y/D = 7.5, the slip lines merge gradually with the jet boundary, making the fourth shock cell difficult to identify. Meanwhile, the organized large scale turbulent vortices are generated periodically and convected downstream with the flow. Panda^[15] observed that a new shock will occur when the turbulent disturbance convects over an existing shock wave in an axisymmetric jet with an NPR of 2.4, and regarded this phenomenon as shock splitting. For the NPR examined in this study, the shocks inside the jets are generally at a quasi-steady state, only oscillating in a limited spatial region without splitting.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) (a) Time-averaged and (b)–(g) instantaneous density gradient at NPR of 5.60 corresponding to t/t0 of (b) 468, (c) 470, (d) 472, (e) 474, (f) 476, and (g) 478, with one oscillation cycle.

Besides the noticeable wavy motion shown in Fig. 7, a careful examination of numerical schlieren photographs (i.e. density gradient) from frame to frame shows a vertical up and down motion for the shocks near the nozzle exit. In particular, figure 8 presents a close-up view of instantaneous flow structures in the near nozzle region at different time for an NPR = 9.34. The vertical oscillation of the first normal shock (i.e., Mach disk) can be identified. The second normal shock oscillates in the same manner as the Mach disk, being up and down alternately. Note that the initial formation and development of the near-nozzle shock structures and Mach disk in highly underexpanded jets have been examined in previous studies, see for example the measurements by Rogers et al.^[24] and the recent LES modelings by Vuorinen et al.^[27] and Hamzehloo and Aleiferis.^[29,30] However, the vertical oscillation of Mach disk when the jet is fully developed, as shown in Fig. 8, is rarely reported. In addition, figures 7 and 8 each indicate that the shock oscillation amplitude is relatively small, implying that it is difficult to measure accurately by, for example, schlieren photographs. On the other hand, more details of the shock motion in underexpanded jets could be achieved based on the high-resolution LES data, especially the quantitative data to be addressed.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Instantaneous contours of density gradient for NPR = 9.34, showing the Mach disk is located at the highest, middle, and lowest position during one oscillation cycle in panels (a)–(c), respectively.

The key to quantitatively investigating the shock oscillation phenomenon is to accurately measure the position of the shocks and the distance over which they move. The shock is generally characterized by an intense jump in pressure as shown in Fig. 6(a), thus the instantaneous profiles of pressure along the jet centerline at different time are examined to quantify the Mach disk oscillation. Based on the detailed analysis of the pressure data, it is found that the Mach disk oscillates at the same cycle around 76t₀ = 190 μs for different NPRs. Figure 9(a) shows the transference process of Mach disk from the lowest to the highest position during half an oscillation cycle for NPR = 5.60. The time-averaged pressure profile is also included in Fig. 9(a), which denotes a much thicker Mach disk than the instantaneous one. In addition, the pressure curves are not uniformly distributed within the oscillation region, but appear to concentrate more near the lowest or highest positions. This indicates that the vertical oscillation of the Mach disk is nonlinear, i.e., moving quickly in the middle of the oscillation region but slowly at the top and bottom end. For other NPRs, the Mach disk oscillates similarly to that with an NPR of 5.60, i.e., undergoing the same oscillation cycle around 76t₀ = 190 μs and moving more quickly in the middle region as shown in Fig. 9(b) for NPR = 9.34.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Instantaneous pressures along the jet centerline, indicating the Mach disk motion in half oscillation cycle at NPR of (a) 5.60 and (b) 9.34.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (color online) Pressure history on either side of the jets at (a) NPR = 5.60, y/D = 2, (b) NPR = 5.60, y/D = 6, (c) NPR = 9.34, y/D = 2, and (d) NPR = 9.34, y/D = 6. The red pressure curves are sampled at x/D = 1 and z/D = 0, while the blue pressure curves are sampled at x/D = −1 and z/D = 0.

The oscillation of Mach disk will bring in pressure fluctuating. Figure 10 shows the time histories of pressure near the jet shear layer at y/D = 2 and 6 for NPR = 5.60 and 9.34, respectively. The pressure fluctuations on either side of jet appear to be correlated and opposite in phase at two different streamwise positions for NPR = 5.60. At a higher NPR of 9.34, the relation between the two pressure signals is complex, opposite or consistent randomly in phase.

Figure 11 shows the spectra and relative phases in the near nozzle region of y/D = 2 for different NPRs. As can be seen, there are several discrete peak frequencies in the pressure spectra. In particular, the frequency peaks circled with solid lines are the shock screech frequencies f_s, and are equal to 37.086 kHz, 34.437 kHz, 31.787 kHz, and 29.801 kHz for NPR = 5.60, 7.47, 9.34, and 11.21, respectively. These screech frequencies decrease with NPR increasing, and have a relative phase angle close to 180° (or −180°), indicating that the dominant modes of jets are all helical.^[35] Besides the peak frequencies of 37.086 kHz and 34.437 kHz, other valid shock screech tones exist in the spectrum for NPR = 5.60 and 7.47. The secondary shock screech tones take values of f_2s = 45.695 kHz and 40.397 kHz for NPR = 5.60 and 7.47, respectively, and are in helical mode as well with a relative phase angel close to 180° (or −180°). In addition, there are two axisymmetric modes with a relative phase angle close to 0° in the spectrum, which takes the same frequency values of f = 5.298 kHz and 14.569 kHz for different NPRs. Remember that the Mach disk oscillation cycles for different NPRs are around T_soc ≈ 76t₀ = 190 μs, which is very close to the cycle of T_f = 1/f = 1/5.298 × 10³ ≈ 188.75 μs for the first axisymmetric mode of 5.298 kHz. Therefore, it can be inferred that the oscillations of Mach disk inside underexpanded jets at different NPRs are dominated by the same axisymmetric mode of 5.298 kHz.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. (color online) Cross spectra and relative phases of pressure fluctuations on either side of jets at NPR of (a) 5.60, (b) 7.47, (c) 9.34, and (d) 11.21.

Note that the multiple screech tones have been observed experimentally in underexpanded rectangular^[18] and circular^[19] jets. In particular, the fully expanded Mach number of the jets with multiple screech tones in Ref. [19] is around 1.8, which is close to those in the NPR = 5.60 and 7.47 cases.^[35] In addition, Raman^[7] pointed out that the symmetric mode produced by a 2D convergent-divergent nozzle may be involved with the nozzle block design. The present LES results presented in the following Subsection 4.4 with different NPRs, grid resolutions, and subgrid scale turbulence models produce the same symmetric modes of f = 5.298 kHz and 14.569 kHz. Therefore, the appearance of the axisymmetric mode of f = 5.298 kHz associated with the Mach disk oscillation at the same frequency for different NPRs is probably related to the current nozzle geometry design.

Figure 12(a) shows the time-averaged as well as instantaneous radical profiles of vorticity magnitude Ω for NPR = 5.60. In particular, the instantaneous curves are consistent with those in Figs. 7(b)–7(g), respectively. As can be seen, the positions at which Ω peaks in the regions of −1.0 < x/D < 0.4 and 0.4 < x/D < 1.0 indicate the jet boundaries. The drifting of peak positions is thus corresponding to the oscillation of shock cells. Therefore, the positions with the maximum value of Ω in the radical direction at different time are carefully identified to quantitatively examine the shock oscillation phenomenon. Figure 13 shows the statistical results for NPR = 5.60.

Fig. 12.

	Figure Option ViewDownloadNew Window
	Fig. 12. (color online) Radial profiles of vorticity magnitude Ω sub-grid scale turbulent kinetic energy ksgs at (a) NPR = 5.60 at y/D = 6, and (b) NPR = 9.34 at y/D = 10. The instantaneous curves in panel (a) are consistent with those in Figs. 7(b)–7(g), successively.

Fig. 13.

	Figure Option ViewDownloadNew Window
	Fig. 13. (color online) Radial positions of maximum value of Ω, revealing shock oscillation phenomenon for NPR = 5.60 at y/D of (a) 6 and (b) 8. The empty circles with black solid lines correspond to those in Figs. 7(b)–7(g) as well as the instantaneous curves in Fig. 12(a), successively.

From Fig. 13, the shock oscillates almost sinusoidally around the time-averaged position, which is in accordance with the antisymmetric screech mode of the jet. Panda^[15] simply used the screech cycle to describe the shock motion versus time. André et al.^[20] demonstrated that the shocks oscillate at the screech frequency by comparing the spectra of the near-field microphone signal with the spectra of the signal containing the shock locations. Figure 13 shows that most of shock oscillation cycles contain six points and the corresponding cycle is about T₁ ≈ 10t₀ = 25.0 μs, which is close to one wave period of the shock screech period 1/f_s = 1/(37.086 × 10³) ≈ 26.96 μs. This demonstrates the direct evidence for the coincidence between the shock oscillation cycle and the screech frequency. In addition, a close examination shows that there are some shock oscillation cycles that contain only five points, implying a shorter cycle period of about T₁ ≈ 8t₀ = 20.0 μs, which is close to one wave period of the secondary shock screech 1/f_2s = 1/(45.695 × 10³) ≈ 21.88 μs. This indicates that the shock cells inside the jet alternately oscillate at two different screech frequencies of f_s and f_2s. The small differences between the oscillation cycle period and the wave period of the screech tone may be caused by insufficient time sampling rate (the instantaneous LES results are saved every 2t₀ = 5.0 μs). The external shock oscillation reflects the internal helical motion of the jet core. From the present finding, it can be inferred that both the peak frequencies of f_s and f_2s are the valid screech frequencies that correspond to different helical modes, and the jet switches the motion modes between the two helical modes randomly for NPR = 5.60 and 7.47 jets, as the conclusion^[35] obtained previously.

Fig. 14.

	Figure Option ViewDownloadNew Window
	Fig. 14. (color online) (a) Amplitude of screech frequencies, and (b) amplitude ratio between the axisymmetric mode of f = 5.298 kHz and screech tone. The screech frequencies presented here are 37.086 kHz, 34.437 kHz, 31.787 kHz, and 29.801 kHz for NPR = 5.60, 7.47, 9.34, and 11.21, respectively.

Figure 13 also shows that the shock oscillation amplitude is about 0.1D at y/D = 6, and grows to around 0.2D at y/D = 8. This quantitatively confirms the previous observation based on Fig. 7 that the oscillation increases for the shocks formed in the further down stream from the nozzle exit. It also agrees reasonably with the measurement by Panda.^[15] Figure 14(a) shows the amplitude of screech frequency, Ψ_s, which increases with the distance of streamwise position increasing, as the shock oscillation magnitude. Meanwhile, the shock oscillation amplitude ∼ 0.2D at y/D = 8 for NPR = 5.60 is larger than ∼ 0.08D for NPR = 9.34 (as shown in Fig. 15), which is consistent with the screech amplitude trend, i.e., the screech amplitude Ψ_s is larger for the lower NPR. These observations confirm the conclusions by André et al.^[20] based on the schlieren photographs, i.e., the shock oscillation magnitude is strongly related to the screech amplitude, probably being an increasing function of the screech amplitude.

Moreover, it is interesting to find from Fig. 13 that the shocks tend to oscillate more intensively toward the jet centerline. First, the shock oscillation amplitude on the centerline side appears to be larger than that away from the jet centerline, no matter whether the value of y/D is 6 or 8, and no matter whether the half shaft is the positive half shaft (x/D > 0) or negative half shaft (x/D < 0). Second, there are generally more time sampling points (usually 3–4) on the centerline side during one oscillation cycle. This finding is in accordance with the observation by Panda,^[15] that is, the shocks inside the underexpanded jets move most in the jet core and least in the shear layer spatially.

Besides the vorticity magnitude Ω, the sub-grid scale turbulent kinetic energy k^sgs also peaks at the jet boundary, and can be used to track the shock locations as well. Figure 12(b) shows the instantaneous radial profiles of k^sgs at different time, and figure 15 presents the jet boundaries identified based on k^sgs for NPR = 9.34. The overall characteristics of shock oscillation for NPR = 9.34 are similar to those for NPR = 5.60. For example, the shock oscillation amplitude is larger for higher streamwise position. It is also found that the shocks appear to oscillate more intensively toward the jet centerline. However, the periodical characteristic of shock oscillation for NPR = 9.34 is more complex than that for NPR = 5.60. There seems to be more concentrated time sampling points near the time-averaged shock positions. The shock also appears to oscillate more quickly for NPR = 9.34 than for NPR = 5.60. These may have two origins. First, the axisymmetric mode plays a more important role and takes effect within a much longer distance at a higher NPR. Figure 14(b) presents the amplitude ratios between axisymmetric mode and screech tone, Ψ_{5.298 kHz}/Ψ_s, for different NPRs. As can be seen, the amplitude ratio Ψ_{5.298 kHz}/Ψ_s is generally large in the near-nozzle region and decreases with being away from the streamwise position. This agrees reasonably with the findings by Mattingly^[50] on low Reynolds jets based on the inviscid stability theory, i.e., the dominant disturbance in the very near field of jets is an axisymmetric one (azimuthal wavenumber m = 0), while the single helical disturbance (|m = 1) is dominant in the further down stream. Meanwhile, the amplitude ratio Ψ_5.298 kHz/Ψ_s is larger for higher NPR at different streamwise positions. These observations suggest that the axisymmetric mode affects the oscillation of shock cells in radial direction except its domination in Mach disk oscillation in the streamwise direction, especially for higher NPR. As a result, the shock oscillation for higher NPR looks more symmetric, which implies that more time sampling points are located near the time-averaged shock positions. Second, an increase of the effect of helical modes with higher frequency and larger wave numbers, such as the double helical modes (|m| = 2), on the jet flow for higher NPR, may result in a faster switching in shock oscillation cycle.

Fig. 15.

	Figure Option ViewDownloadNew Window
	Fig. 15. (color online). Radial positions of maximum value of ksgs, revealing the shock oscillation phenomenon for NPR = 9.34 at y/D of (a) 8 and (b) 10. The empty circles with black solid lines correspond to the instantaneous curves in Fig. 12(b).

In general, the fine mesh is able to capture more small-scale turbulence structures, but its computational cost is higher. Likewise, the one-equation SGS model has less advantage in computational efficiency because more equations are solved than the zero equation models, such as the Smagorinsky model. Therefore, the Smagorinsky model^[51,52] is also used for simulating the NPR = 5.60 jet, and the effect of grid resolution and SGS model on revealing the unsteady shock motion are investigated in this section.

Fig. 16.

	Figure Option ViewDownloadNew Window
	Fig. 16. (color online) Instantaneous pressures along the jet centerline obtained under different simulation conditions for NPR = 5.60, indicating the Mach disk motion in half the oscillation cycle, by using (a) one equation SGS model with coarse mesh and (b) Smagorinsky model with fine mesh.

Figures 16(a) and 16(b) show the Mach disk oscillations for NPR = 5.60 jet, revealed using the coarse mesh and the Smagorinsky model, respectively, where the time is the same as that presented in Fig. 9(a). Meanwhile, like Fig. 11(a), the cross spectrum and relative phase of pressure fluctuation at y/D = 2 obtained by the coarse mesh and the Smagorinsky model are plotted in Fig. 17. As can be seen, the oscillation characteristics of the Mach disk are also well captured by the coarse mesh and the Smagorinsky model, i.e., the Mach disk oscillates at the same cycle around 76t₀ = 190 μs and moves more quickly in the middle region. Moreover, the axisymmetric mode of 5.298 kHz that dominates the Mach disk oscillation also exists in the spectra obtained by the coarse mesh and Smagorinsky model, which may be related to the nozzle geometry used in the present LES as discussed previously. However, the oscillation amplitude of Mach disk presented in Fig. 16 seems to be smaller than that shown in Fig. 9(a), especially for the results obtained using the coarse mesh.

Fig. 17.

	Figure Option ViewDownloadNew Window
	Fig. 17. (color online) Cross spectra and relative phases of pressure fluctuation on either side of the jets obtained under different simulation conditions for NPR = 5.60 by using (a) one equation SGS model with coarse mesh and (b) Smagorinsky model with fine mesh.

As also shown in Fig. 17, the two helical modes that dominate the shock oscillation in the radial direction are captured by the coarse mesh, and take values of 35.761 kHz and 44.370 kHz, respectively, which are slightly smaller than 37.086 kHz and 45.695 kHz resolved by the fine mesh. Meanwhile, the Smagorinsky model captures the two helical modes as well, and the peak frequencies are the same as those obtained by the one-equation model, i.e., 37.086 kHz and 45.695 kHz.

Fig. 18.

	Figure Option ViewDownloadNew Window
	Fig. 18. (color online) Comparison of shock screech amplitude obtained under different simulation conditions for NPR = 5.60: (a) amplitude of shock screech frequency f2s, (b) non-dimensional amplitude ratio of shock screech tone f2s.

In the previous section, it is indicated that the magnitude of shock oscillation in the radial direction is an increasing function of the screech amplitude. Therefore, the amplitudes of shock screech frequency f_2s signals captured using different simulation conditions are compared in Fig. 18(a) to characterize the radial shock oscillation magnitude. As can be seen, the overall profiles of the shock screech amplitudes obtained based on the coarse mesh and Smagorinsky model are similar to those in the baseline case (the simulation with the fine mesh and one-equation model) as well as the results presented in Fig. 14 (a), i.e., the screech amplitude Ψ_2s increases with the distance of streamwise position decreasing. However, the simulations using coarse mesh and Smagorinsky model produce relatively small Ψ_2s values. In particular, the Ψ_2s values for the coarse mesh take the smallest values at y/D = 6 and 8. This observation suggests that the coarse mesh tends to predict a much weaker radial shock oscillation, just as predicting a weaker Mach disk oscillation shown in Fig. 16. In order to further check how far it deviates from the screech amplitude obtained based on the one-equation model with fine mesh, the screech amplitude of the baseline case is chosen as the reference value, and the non-dimensional amplitude ratio Ψ_2s/Ψ_baseline using the Smagorinsky model with coarse mesh is computed numerically, and the results are plotted in Fig. 18(b). It is clear that the amplitude ratio for the coarse mesh decreases with the distance of streamwise position increasing. This trend may be explained by the increase of dissipation due to the decay of grid resolution along the streamwise direction. For the simulation using the Smagorinsky model, the amplitude ratio is about 0.68 at y/D = 2, and increases gradually to 0.74 at y/D = 8. These findings indicate that the disadvantage of Smagorinsky model in revealing the unsteady shock motions in underexpanded jets is not so remarkable, because of the considerable performances in capturing the frequency and magnitude of shock oscillation. On the other hand, it appears that the shock oscillation characteristics resolved by the LES, including the oscillation frequency and magnitude, rely more on the grid resolution used.